
#include "multimin/gsl_ext_function.h"
#include "multimin/gsl_ext_multimin.h"
#include "multimin/gsl_ext_testfunctions.h"

#include <gsl/gsl_multimin.h>
#include <stdlib.h>

static const char *ROSENBROCK_EXPR = "100*(y-x^2)^2+(1-x)^2";

static void print_results(const char *algorithm_name,
                          int n_iter,
                          const gsl_vector *x,
                          double f,
                          int status)
{
  if(status == GSL_SUCCESS)
    printf("%-15s: niter: %-4d x*=%.5f  y*=%.5f  f(x*,y*)=%.5e\n",
           algorithm_name,
           n_iter,
           gsl_vector_get(x, 0),
           gsl_vector_get(x, 1),
           f);
  else
    printf("%-15s: FAILURE (%s)\n", algorithm_name, gsl_strerror(status));
}

static void test_fminimizers(gsl_ext_func_eval_type f_eval_type)
{
  const gsl_multimin_fminimizer_type *T = 
    gsl_ext_multimin_fminimizer_lrwwsimplex;
  gsl_multimin_fminimizer *s;
  gsl_multimin_function *func;
  size_t iter = 0;
  int status;
  gsl_vector *x0;
  
  x0 = gsl_vector_alloc(2);
  gsl_vector_set(x0, 0, -1.2);
  gsl_vector_set(x0, 1, 1.0);
  
  gsl_ext_testfunction_params f_params = { 2, 2 };
  
  if(f_eval_type == GSL_EXT_FUN_EVAL_SYM)
    func = gsl_ext_multimin_function_alloc(ROSENBROCK_EXPR,
                                           GSL_EXT_FUN_EVAL_SYM,
                                           &f_params);
  else
    func = gsl_ext_multimin_function_alloc("extendedrosenbrock",
                                           GSL_EXT_FUN_EVAL_PRECOMP,
                                           &f_params);
  s = gsl_multimin_fminimizer_alloc(T, 2);
  
  gsl_vector *iv = gsl_vector_alloc(2);
  gsl_vector_set(iv, 0, 1.0);
  gsl_vector_set(iv, 1, 1.0);
  gsl_multimin_fminimizer_set(s, func, x0, iv);
  gsl_vector_free(iv);
  
  do
  {
    iter++;
    status = gsl_multimin_fminimizer_iterate(s);
    
    if(status)
      break;
    
    status = gsl_multimin_test_size(s->size, 1e-14);
  }
  while(status == GSL_CONTINUE && iter < 10000);
  
  print_results(T->name, iter, s->x, s->fval, status);
  
  gsl_ext_multimin_function_free(func);
  gsl_multimin_fminimizer_free(s);
  gsl_vector_free(x0);
}

static void test_fdfminimizers(gsl_ext_func_eval_type f_eval_type)
{
  const gsl_multimin_fdfminimizer_type *ALGORITHM_TYPES[] = 
  {
    gsl_ext_multimin_fdfminimizer_bfgs_hess_f,
    gsl_ext_multimin_fdfminimizer_bfgs_f,
    gsl_ext_multimin_fdfminimizer_bfgs_hess_mt,
    gsl_ext_multimin_fdfminimizer_bfgs_mt,
    gsl_ext_multimin_fdfminimizer_conjgrad_fr_mt,
    gsl_ext_multimin_fdfminimizer_conjgrad_pr_mt
  };
  const int NUM_ALGORITHMS = 6;
  
  gsl_multimin_fdfminimizer *s;
  gsl_multimin_function_fdf *f;
  int i;
  size_t iter = 0;
  int status;
  gsl_vector *x0;
  
  x0 = gsl_vector_alloc(2);
  gsl_vector_set(x0, 0, -1.2);
  gsl_vector_set(x0, 1, 1.0);
  
  gsl_ext_testfunction_params f_params = { 2, 2 };
  
  if(f_eval_type == GSL_EXT_FUN_EVAL_SYM)
    f = gsl_ext_multimin_function_fdf_alloc(ROSENBROCK_EXPR,
                                            GSL_EXT_FUN_EVAL_SYM,
                                            GSL_EXT_DERIV_EVAL_SYM,
                                            &f_params);
  else
    f = gsl_ext_multimin_function_fdf_alloc("extendedrosenbrock",
                                            GSL_EXT_FUN_EVAL_PRECOMP,
                                            GSL_EXT_DERIV_EVAL_FD,
                                            &f_params);
  
  for(i = 0; i < NUM_ALGORITHMS; i++)
  {
    s = gsl_multimin_fdfminimizer_alloc(ALGORITHM_TYPES[i], 2);
    gsl_multimin_fdfminimizer_set(s, f, x0, 1.0, 0.1);
    
    iter = 0;
    do
    {
      iter++;
      status = gsl_multimin_fdfminimizer_iterate(s);
      
      if(status)
        break;
      
      status = gsl_multimin_test_gradient(s->gradient, 1e-8);
    }
    while(status == GSL_CONTINUE && iter < 1000);
    
    print_results(ALGORITHM_TYPES[i]->name, iter, s->x, s->f, status);
    gsl_multimin_fdfminimizer_free(s);
  }

  gsl_ext_multimin_function_fdf_free(f);
  gsl_vector_free(x0);
}

static void test_fd2fminimizers(gsl_ext_func_eval_type f_eval_type)
{
  const gsl_ext_multimin_fd2fminimizer_type *T = 
    gsl_ext_multimin_fd2fminimizer_mnewton;
  gsl_ext_multimin_fd2fminimizer *s;
  gsl_ext_multimin_function_fd2f *func;
  
  size_t iter = 0;
  int status;
  
  gsl_vector *x0;
  
  x0 = gsl_vector_alloc(2);
  gsl_vector_set(x0, 0, -1.2);
  gsl_vector_set(x0, 1, 1.0);
  
  gsl_ext_testfunction_params func_params = { 2, 2 };
  if(f_eval_type == GSL_EXT_FUN_EVAL_SYM)
    func = gsl_ext_multimin_function_fd2f_alloc(ROSENBROCK_EXPR,
                                                GSL_EXT_FUN_EVAL_SYM,
                                                GSL_EXT_DERIV_EVAL_SYM,
                                                GSL_EXT_DERIV_EVAL_SYM,
                                                &func_params);
  else
    func = gsl_ext_multimin_function_fd2f_alloc("extendedrosenbrock",
                                                GSL_EXT_FUN_EVAL_PRECOMP,
                                                GSL_EXT_DERIV_EVAL_FD,
                                                GSL_EXT_DERIV_EVAL_FD,
                                                &func_params);
  s = gsl_ext_multimin_fd2fminimizer_alloc(T, 2);
  gsl_ext_multimin_fd2fminimizer_set(s, func, x0, 1.0, 0.1);
  
  do
  {
    iter++;
    status = gsl_ext_multimin_fd2fminimizer_iterate(s);
    
    if(status)
      break;
    
    status = gsl_multimin_test_gradient(s->gradient, 1e-6);
  }
  while(status == GSL_CONTINUE && iter < 1000);
  
  print_results(T->name, iter, s->x, s->f, status);
  
  gsl_ext_multimin_function_fd2f_free(func);
  gsl_ext_multimin_fd2fminimizer_free(s);
  gsl_vector_free(x0);
}

int main(int argc, char *argv[])
{
  printf("Testing GSL++ algorithms on the extended Rosenbrock function (symbolic):\n");
  test_fminimizers(GSL_EXT_FUN_EVAL_SYM);
  test_fdfminimizers(GSL_EXT_FUN_EVAL_SYM);
  test_fd2fminimizers(GSL_EXT_FUN_EVAL_SYM);
  printf("Testing GSL++ algorithms on the extended Rosenbrock function (precompiled):\n");
  test_fminimizers(GSL_EXT_FUN_EVAL_PRECOMP);
  test_fdfminimizers(GSL_EXT_FUN_EVAL_PRECOMP);
  test_fd2fminimizers(GSL_EXT_FUN_EVAL_PRECOMP);
  
  return EXIT_SUCCESS;
}
